Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixes and improvements in SOLPSSimulation, SOLPSMesh and related functions #35

Merged

Conversation

vsnever
Copy link
Member

@vsnever vsnever commented Sep 4, 2020

This PR addresses the issues #8, #20, #31, #33, #34, #36 and #37.

Issue #8 is addressed in the following way:

  • Safe to use SOLPSFunction2D and SOLPSVectorFunction2D classes are implemented. They replace SOLPSFunction3D and SOLPSVectorFunction3D. For 3D, use AxisymmetricMapper(SOLPSFunction2D) and VectorAxisymmetricMapper(SOLPSVectorFunction2D).

  • SOLPSMesh and SOLPSSimulation can now be initialised without modification of hidden attributes;

  • SOLPSMesh has new attributes for basis vectors, cell connection areas, indices of neighbouring cells and new functions to_poloidal() and to_cartesian() for converting vectors defined on a grid from/to (poloidal, radial)/(R, Z);

  • species_list in SOLPSSimulation is a list of (name, charge) pairs now instead of a list of strings;

  • Setting b_field and velocities vectors in curvilinear poloidal coordinates (epol, erad, etor) automaticaly sets b_field_cylindrical and velocities_cylindrical (eR, eφ, eZ) and vice versa (note: etor = -eφ);

  • Incorrect calculation of basis vectors using cell centres is fixed (basis vectors are defined by left and bottom faces of the cell);

  • The code for mesh.triangles and mesh.triangle_to_grid_map creation is vectorised (as all other code wherever possible);

  • SOLPSSimulation now has setters for its properties;

  • load_solps_from_mdsplus(), load_solps_from_raw_output() and load_solps_from_balance() do not brake SOLPSSimulation interface anymore;

  • load_solps_from_mdsplus() and load_solps_from_raw_output() are made more robust. Now both support B2 standalone runs. Error handling is added to load_solps_from_mdsplus() for some optional attributes that may not be present on the server;

  • Other small improvements.

Issue #20 is fixed. *_f2d() and *_f3d() interpolators are added for plasma parameters. For vectors _carteisan suffix is used instead of _f3d.

Issue #31 is fixed. Now if the neutral densities from Eirene are available, they replace the neutral densities from B2.

Issue #33 is fixed using lookup_isotope() and lookup_element() methods. The plasma composition is defined from nuclear charge, atomic mass number and electrical charge. This is more robust than parsing a string with species names.

Issue #34 is addressed by parsing the following additional quantities in load_solps_from_mdsplus() and load_solps_from_raw_output():

  • Ion temperature;
  • magnetic field (previously was read only from MDS);
  • atom/ion fluxes (B2), neutral atom densities (Eirene, if available),
  • neutral atom fluxes (Eirene, if available);
  • neutral atom temperatures (Eirene, if available).

The fluxes are not stored in SOLPSSimulation, but used to calculate velocities.

Issue #36 is fixed. The parallel velocity and the B2 fluxes (in s-1) defined on cell faces are used in b2_flux_to_velocity() to calculate velocities at cell centres. If neutral atom fluxes from Eirene (in m-2 s-1) defined at cell centres are available, they are used in eirene_flux_to_velocity() replacing the B2 ones.

Issue #37 is fixed. Now all arrays are row-major. However, the indexing is inverted. The correct index order now is: (radial, poloidal) or (species, radial, poloidal) for scalars and (vector component, radial, poloidal) or (species, vector component, radial, poloidal) for vectors. This should improve the performance in typical workflows.

Remaining issues:

  • Electron velocities are still not read from SOLPS data.

  • Molecular and total H-alpha emission are still not read from SOLPS data.

@vsnever
Copy link
Member Author

vsnever commented Sep 5, 2020

Just updated the code by adding ion and neutral atom temperatures. The latter is not available in stand-alone B2.5 simulations.

@vsnever vsnever changed the title Add setters to SOLPSSimulation, remove unsued attributes, improved reading from MDSPlus and raw files Fixes and improvements in SOLPSSimulation, SOLPSMesh and related functions Sep 5, 2020
@vsnever
Copy link
Member Author

vsnever commented Sep 5, 2020

I renamed this PR because now it addresses not just #31 but in some degree #8 too.

Replacing the algorithm for poloidal basis vectors calculation in SOLPSMesh was more a bug fix than an improvement. The SOLPS mesh have break points due to the separatrix. In the old algorithm, which calculated the basis vectors using the cell centres, the basis vectors had wrong values at these break points.

@vsnever
Copy link
Member Author

vsnever commented Sep 6, 2020

Just noticed that conversion from poloidal to Cartesian coordinates in mdsplus.py is incorrect. It's for orthogonal basis, but SOLPS uses curvilinear coordinates.

@vsnever
Copy link
Member Author

vsnever commented Sep 6, 2020

Just noticed that conversion from poloidal to Cartesian coordinates in mdsplus.py is incorrect. It's for orthogonal basis, but SOLPS uses curvilinear coordinates.

Sorry, I'm wrong here. It converts to Cartesian, so it's OK.

@mattngc
Copy link
Member

mattngc commented Sep 6, 2020

Hey Vlad, I don't have any strong opinions on these improvements. I don't use this module as much at the moment, so I think opinions from the main users are more important. I'll let @jacklovell and @Mateasek give feedback on these changes.

…_simulation_files. Removed radial_particle_flux and radial_area attributes from SOLPSSimulation.
@vsnever
Copy link
Member Author

vsnever commented Sep 7, 2020

Hey Vlad, I don't have any strong opinions on these improvements. I don't use this module as much at the moment, so I think opinions from the main users are more important. I'll let @jacklovell and @Mateasek give feedback on these changes.

Hi Matthew, thanks for reply. I need a couple more days to get this PR ready for review.

@vsnever vsnever marked this pull request as ready for review September 9, 2020 13:28
@vsnever
Copy link
Member Author

vsnever commented Sep 9, 2020

In the first post, I summaraised the changes made in this PR and marked it ready for review.

Hi @jacklovell and @Mateasek, I am awaiting your verdict).

Hello @jrh-ccfe, could you please share any balance.nc file, so I could improve load_solps_from_balance(). I have access to Heimdall cluster, so a path to the file will be enough.

@vsnever
Copy link
Member Author

vsnever commented Sep 15, 2020

I updated balance.py and revert the change in the fort44 parser.

Copy link
Member

@jacklovell jacklovell left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this Vlad: it's a huge amount of great work and a big improvement to the module overall. I've made some in-line comments on the code, but generally it is of a very good standard.

I've run the 3 demo scripts we currently have on a limited set of runs that I have access to and they seem to produce mostly the same results as before. I expect any small differences are a result of bug fixes. Have you done regression testing on any other SOLPS cases you have?

There's a few places in the code with large commented-out blocks: these can be removed if they won't be added back in later. I'm thinking in particular in solps_plasma.py, where there are several commented-out plotting methods which should live in demo scripts instead.

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

cherab/solps/b2/parse_b2_block_file.py Outdated Show resolved Hide resolved
cherab/solps/b2/parse_b2_block_file.py Outdated Show resolved Hide resolved
cherab/solps/eirene/parser/fort44_2013.py Outdated Show resolved Hide resolved
cherab/solps/eirene/parser/fort44_2017.py Outdated Show resolved Hide resolved
cherab/solps/eirene/eirene.py Show resolved Hide resolved
if self._total_rad is None:
raise RuntimeError("Total radiation not available for this simulation")
if self._total_radiation is None:
raise RuntimeError("Total radiation not available for this simulation.")
else:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

else is not necessary after a raise. This applies to all the total_radiation*, b_field* and eirene_simulation properties.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code was before me, some properties threw errors when not initialized and some did not. I think, they should not throw errors. If a property is required for a function to work, that function must check if the value is set and throw a RuntimeError if not. I've removed these raise statements from properties and instead added them to functions that don't work without these properties.

cherab/solps/solps_plasma.py Outdated Show resolved Hide resolved
for k, sp in enumerate(self.species_list):
if self.electron_velocities_cartesian is None:
print('Warning! No electron velocity data available for this simulation.')
electron_velocity = lambda x, y, z: Vector3D(0, 0, 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be a ConstantVector3D object, not a lambda function.

if isinstance(self.velocities_cartesian, np.ndarray):
velocity = SOLPSVectorFunction3D(tri_index_lookup, tri_to_grid, self.velocities_cartesian[:, :, k, :])
if self.velocities_cartesian is not None:
velocity = self.velocities_cartesian[k]
else:
velocity = lambda x, y, z: Vector3D(0, 0, 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be a ConstantVector3D object, not a lambda function.

else:
electron_velocity = self.electron_velocities_cartesian
plasma.electron_distribution = Maxwellian(self.electron_density_f3d, self.electron_temperature_f3d, electron_velocity, electron_mass)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Trailing whitespace.

@jacklovell
Copy link
Member

One other comment: if Axisymmetric wrappers of SOLPSFunction2D and SOLPSVectorFunction2D replace the corresponding 3D functions, perhaps we should remove the latter from the module entirely. The 3D functions are still present in this PR.

@CnlPepper
Copy link
Member

One other comment: if Axisymmetric wrappers of SOLPSFunction2D and SOLPSVectorFunction2D replace the corresponding 3D functions, perhaps we should remove the latter from the module entirely. The 3D functions are still present in this PR.

I agree, it would be best to remove redundant code.

@CnlPepper
Copy link
Member

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

While there are large API changes to the SOLPSSimulation API, the core functionality for loading and generating Plasma objects
and the object layout are essentially the same. It is definitely a new minor release, possibly a major release.

The arguments for a major/minor version increment would be the degree of impact to the main usage of the package. This one lies on the border I suspect. As the active users of this package, you would be best placed to assess the impact on users and then decide if you want to go for v2.0.0 or v1.2.0.

Please do add a changelog. I would recommend getting into the habit of checking if it needs updating with every pull request.

must be 3 dimensional. Example shape is (4 x 32 x 98).
In SOLPS notation: left/right - poloidal prev./next, bottom/top - radial prev./next.
Cell indexing starts with 0 and -1 means no neighbour.
:param ndarray neighbix: Array of radial indeces of neighbouring cells in order: left, bottom, right, top,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't it be neighbiy?


self._eirene = value

def b2_flux_to_velocity(self, poloidal_flux, radial_flux, parallel_velocity):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't it be better to make this method a standalone instead of class method, like it is in the case of EFITEquilibrium class? This class is becoming quite big, it would help decreasing it. Also the method itself could be useful/used for various purposes. The same could be done with eirene_flux_to_velocity.

@Mateasek
Copy link
Member

Mateasek commented Sep 18, 2020

That is a large amount of work @vsnever, very nice. Thanks for all the work you put into it, especially when SOLPS is not the easiest code to understand and read! I would like to propose one more change, since this is quite major rewrite. I would suggest also removing the commented plotting methods from the SOLPSimulation class. I think the plots would be better done outside in a separate file, as it is with equilibrium plots in Cherab. I don't mean that you should be also adding it now.

…inside/outside mesh in SOLPSTotalRadiatedPower.
@vsnever
Copy link
Member Author

vsnever commented Sep 19, 2020

@jacklovell, @Mateasek, @CnlPepper, thank you very much for the code review. I hope I addressed all the issues you raised. For three of them I left the comments (see above). If I did not leave a comment, it means that the issue is fixed exactly as you suggested.

I've run the 3 demo scripts we currently have on a limited set of runs that I have access to and they seem to produce mostly the same results as before. I expect any small differences are a result of bug fixes. Have you done regression testing on any other SOLPS cases you have?

When the data is read from the MDSPlus, the only difference should be in species velocities, because of the bug fix. When the data is read from the balance.nc file, then in addition to velocities, the density of impurity neutral atoms is different, because now it's obtained from the Eirene data. If the data is read from raw simulation files then the density of main plasma neutral atoms is also different because previously it was obtained from B2.

Also, considering the amount of change here it would be worth adding to the changelog while the changes are still fresh. I'll defer to @CnlPepper on whether this is a big enough change to warrant a 1.2 (or even a 2.0) release, or whether a 1.1.1 would be suitable.

I updated the changelog suggesting a 1.2 version.

In addition to the improvements you suggested, I removed redundant check for inside/outside mesh in SOLPSTotalRadiatedPower emitter, and because now this emitter uses only a total_radiation_f3d from SOLPSSimulation, I think that passing the simulation object to this emitter is no longer needed, a Function3D instance is enough.

Copy link
Member

@jacklovell jacklovell left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me.

@jacklovell
Copy link
Member

@CnlPepper @Mateasek are you happy for these changes to be merged?

@Mateasek
Copy link
Member

Yes, sorry I am new to this, forgot to click on approve. Very nice piece of work!

@CnlPepper
Copy link
Member

Happy here, nice work!

@vsnever
Copy link
Member Author

vsnever commented Sep 21, 2020

@jacklovell, @Mateasek, @CnlPepper, thanks a lot!

@jacklovell, please do not close #34 when this is merged. There is still plenty of work to do.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
5 participants